Critical behavior of magnetic thin films as a function of thickness 
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We study the critical behavior of magnetic thin hims as a function of the him thickness. We 
use the ferromagnetic Ising model with the high-resolution multiple histogram Monte Carlo (MC) 
simulation. We show that though the 2D behavior remains dominant at small thicknesses, there is a 
systematic continuous deviation of the critical exponents from their 2D values. We observe that in 
the same range of varying thickness the deviation of the exponent v is rather small, while exponent 
(3 suffers a larger deviation. We explain these deviations using the concept of "effective" exponents 
suggested by Capehart and Fisher in a finite-size analysis. The shift of the critical temperature with 
the film thickness obtained here by MC simulation is in an excellent agreement with their prediction. 
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I. INTRODUCTION 

During the last 30 years, physics of surfaces and ob- 
jects of nanometric size have attracted an immense inter- 
est. This is due to important applications in industryji^ 
An example is the so-called giant magneto-resistance 
(GMR) used in data storage devices, magnetic sensors, 
..i ?i 4 i?it In parallel to these experimental developments, 
much theoretical effort 7 -^ has also been devoted to the 
search of physical mechanisms lying behind new proper- 
ties found in nanoscale objects such as ultrathin films, 
ultrafme particles, quantum dots, spintronic devices etc. 
This effort aimed not only at providing explanations for 
experimental observations but also at predicting new ef- 
fects for future experiments. 

The physics of two-dimensional (2D) systems is very 
exciting. Some of those 2D systems can be exactly solved: 
one famous example is the Ising model on the square 
lattice which has been solved by Onsager^ This model 
shows a phase transition at a finite temperature T c given 
by sinh 2 (2J/fcsT c ) = 1 where J is the nearest-neighbor 
(NN) interaction. Another interesting result is the ab- 
sence of long-range ordering at finite temperatures for 
the continuous spin models (XY and Heisenberg models) 
in 2Dfi^ In general, three-dimensional (3D) systems for 
any spin models cannot be unfortunately solved. How- 
ever, several methods in the theory of phase transitions 
and critical phenomena can be used to calculate the crit- 
ical behaviors of these systems^ 

This paper deals with systems between 2D and 
3D. Many theoretical studies have been devoted to 
thermodynamic properties of thin films, magnetic 
multilayers, ..i 7 - ^ 12 ' 13 ! 14 In spite of this, several points are 
still not yet understood. It is known a long time ago 
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that the presence of a surface in magnetic materials can 
give rise to surface spin-waves which are localized in the 
vicinity of the surfaced These localized modes may be 
acoustic with a low-lying energy or optical with a high en- 
ergy, in the spin-wave spectrum. Low-lying energy modes 
contribute to reduce in general surface magnetization at 
finite temperatures. One of the consequences is the sur- 
face disordering which may occur at a temperature lower 
than that for interior magnetization^ The existence of 
low-lying surface modes depends on the lattice structure, 
the surface orientation, the surface parameters, surface 
conditions (impurities, roughness, ...) etc. There are two 
interesting cases: in the first case a surface transition oc- 
curs at a temperature distinct from that of the interior 
spins and in the second case the surface transition co- 
incides with the interior one, i. e. existence of a single 
transition. Theory of critical phenomena at surfaces 7,8 
and Monte Carlo (MC) simulation s 17 ' 18 of critical be- 
havior of the surface-layer magnetization at the extraor- 
dinary transition in the three-dimensional Ising model 
have been carried out. These works suggested several 
scenarios in which the nature of the surface transition 
and the transition in thin films depends on many factors 
in particular on the symmetry of the Hamiltonian and on 
surface parameters. 

We confine ourselves here in the case of a simple cu- 
bic film with Ising model. For our purpose, we suppose 
all interactions are the same everywhere even at the sur- 
face. This case is the simplest case where there is no 
surface-localized spin-wave modes and there is only a 
single phase transition at a temperature for the whole 
system (no separate surface phase transition)j 15 i 16 Other 
complicated cases will be left for future investigations. 
However, some preliminary discussions on this point for 
complicated surfaces have been reported in some of our 
previous papers * 19 ' 20 In the case of a simple cubic film 
with Ising model, Capehart and Fisher have studied the 
critical behavior of the susceptibility using a finite-size 
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scaling analysis.™ They showed that there is a crossover 
from 2D to 3D behavior as the film thickness increases. 
The so-called "effective" exponent 7 has been shown to 
vary according to a scaling function depending both on 
the film thickness and the distance to the transition tem- 
perature. As will be seen below the scaling suggested by 
Capehart and Fisher is in agreement with what we find 
here using extensive MC simulation. 

The aim of this paper is to investigate the effect of the 
thickness on the critical exponents of the film. To carry 
out these purposes, we shall use MC simulations with 
highly accurate multiple histogram technique i 22 ' 23 i 24 

The paper is organized as follows. Section II is devoted 
to a description of the model and method. Results are 
shown and discussed in section III. Concluding remarks 
are given in section IV. 



II. MODEL AND METHOD 
A. Model 

Let us consider the Ising spin model on a film made 
from a ferromagnetic simple cubic lattice. The size of 
the film is L x L x N z . We apply the periodic boundary 
conditions (PBC) in the xy planes to simulate an infinite 
xy dimension. The z direction is limited by the film 
thickness N z . If N z = 1 then one has a 2D square lattice. 

The Hamiltonian is given by 



H 



{hi) 



(1) 



where <7, is the Ising spin of magnitude 1 occupying the 
lattice site i, Y^Uj) indicates the sum over the NN spin 
pairs Oi and <jj. 

In the following, the interaction between two NN sur- 
face spins is denoted by J s , while all other interactions 
are supposed to be ferromagnetic and all equal to J = 1 
for simplicity. Let us note in passing that in the semi- 
infinite crystal the surface phase transition occurs at the 
bulk transition temperature when J s ~ 1.52 J. This point 
is called "extraordinary phase transition" which is char- 
acterized by some particular critical exponentSi 17 i 18 In 
the case of thin films, i. e. N z is finite, it has been the- 
oretically shown that when J s = 1 the bulk behavior is 
observed when the thickness becomes larger than a few 
dozens of atomic layers^ surface effects are insignificant 
on thermodynamic properties such as the value of the 
critical temperature, the mean value of magnetization at 
a given T, ... When J s is smaller than J, surface mag- 
netization is destroyed at a temperature lower than that 
for bulk spinsJ^ The criticality of a film with uniform in- 
teraction, i.e. J s = J, has been studied by Capehart and 
Fisher as a function of the film thickness using a scaling 
analysis 2 ^ and by MC simulationsj 25 ' 26 The results by 
Capehart and Fisher indicated that as long as the film 
thickness is finite the phase transition is strictly that of 



the 2D Ising universality class. However, they showed 
that at a temperature away from the transition temper- 
ature T C (N Z ), the system can behave as a 3D one when 
the spin-spin correlation length £(T) is much smaller than 
the film thickness, i. e. £(T)/N Z < 1. As T gets very 
close to T C (N Z ), £{T)/N Z — > 1, the system undergoes a 
crossover to 2D criticality. We will return to this work 
for comparison with our results shown below. 



B. Multiple histogram technique 

The multiple histogram technique is known to repro- 
duce with very high accuracy the critical exponents of 
second order phase transitions j 22 i 23 ' 24 

The overall probability distribution— at temperature 
T obtained from n independent simulations, each with 
Nj configurations, is given by 



P(E,T) 



Y2=iN j exp[E/k B T J -f j ] 



where 



exp[/ < ]=5^P(f;,T < ) 



(2) 



(3) 



The thermal average of a physical quantity A is then 
calculated by 



in which 



(A(T)) = ^AP(E,T)/z(T), 



z(T)=Y,P(E,T). 



(4) 



(5) 



Thermal averages of physical quantities are thus calcu- 
lated as continuous functions of T, now the results should 
be valid over a much wider range of temperature than for 
any single histogram. 

In MC simulations, one calculates the averaged order 
parameter (M) (M: magnetization of the system), av- 
eraged total energy (E), specific heat C v , susceptibility 
X, first order cumulant of the energy Cjj , and n th order 
cumulant of the order parameter V n for n — 1 and 2. 
These quantities are defined as 
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Let us discuss the case where all dimensions can go to 
infinity. For example, consider a system of size L d where 
d is the space dimension. For a finite L, the pseudo 
"transition" temperatures can be identified by the max- 
ima of C v and Xi ■■■■ These maxima do not in general 
take place at the same temperature. Only at infinite L 
that the pseudo "transition" temperatures of these re- 
spective quantities coincide at the real transition tem- 
perature T c (oo). So when we work at the maxima of 
V n , C' v and Xi we are m f ac t working at temperatures 
away from T c (oo). Let us define the reduced tempera- 
ture which measures the "distance" from T c (oo) by 



t = 



T-T c (oo) 
T c (oo) 



(11) 



This distance tends to zero when all dimensions go to 
infinity. For large values of L, the following scaling rela- 
tions are expected (see details in Ref. \24 



V r r x «I 1/ ", y 2 max cx l 1 /", 
CT X = C + dL a ' v 



and 



X max oc 



£7A 



(12) 
(13) 

(14) 



at their respective 'transition' temperatures T C (L), and 



Cu = Cu[T c {oo))+AL 



-a/u 



C. The case of films with finite thickness 

In the case of a thin film of size L x L x N z , Cape- 
hart and Fisher— have showed that as long as the film 
thickness N z is not allowed to go to infinity, there is a 
2D-3D crossover if one does not work at the real tran- 
sition temperature T C (L = oo,N z ). Following Capehart 
and Fisher, let us define 



t 



T - T C (L = oo, N z ) 



T C (3D) 



(20) 
(21) 



(15) 



where v^d is the 3D v exponent and T C (3D) the 3D crit- 
ical temperature. When x is larger than a value xq, i. 
c. at a temperature away from T C (L — oo,N z ), the sys- 
tem behaves as a 3D one. While when x < Xq, it should 
behave as a 2D one. This crossover was argued from a 
comparison of the correlation length in the z direction to 
the film thickness. As a consequence, if we work exactly 
at T C (L = oo,N z ) we should observe the 2D critical ex- 
ponents for finite N z . Otherwise, we should observe the 
so-called "effective critical exponents" whose values are 
found between those of 2D and 3D cases. This point is 
fundamentally very important. There have been some at- 
tempts to verify it by MC simulations^ but these results 
were not convincing due to their poor MC quality. In the 
following we show with high-precision MC technique that 
the prediction of Capehart and Fisher is really verified. 



T c (oo) OC 



-P/v 



and 



T C {L) = T c {oo) + C A L 



-l/u 



(16) 



(17) 



where A, Co, C\ and Ca are constants. We estimate 
v independently from V' 1 max and V^ max . With this value 
we calculate 7 from ^ max and a from C™ ax . Note that 
we can estimate T c (oo) using the last expression. Then, 
using T c (oo), we can calculate (3 from My/00) • The Rush- 
brooke scaling law a + 2(3 + 7 = 2 is then in principle 
verified. 

Let us emphasize that the expressions Eqs. (|T2l - (|r7|) 
are valid for large L. To be sure that L are large enough, 
one has to allow for corrections to scaling of the form, for 
example, 



x max = B 1 L~ r/u (l + B 2 L- UJ ) 



V 1 



D 1 L 1 / V (1 + D 2 L- U ) 



(18) 
(19) 



where B±, B%, D\ and D 2 are constants and a; is a cor- 
rection exponent. 27 Similar forms exist also for the other 
exponents. Usually, these corrections are extremely small 
if L is large enough as is the case with today's large- 
memory computers. So, in general they do not therefore 
alter the results using Eqs. (|12p ~ (TTT)) . 



III. RESULTS 

The xy linear sizes L — 20, 24, 30, ...,80 have been used 
in our simulations. For N z = 3, sizes up to 160 have been 
used to evaluate corrections to scaling. 

In practice, we use first the standard MC simula- 
tions to localize for each size the transition tempera- 
tures T t f(L) for specific heat and T m (L) for suscepti- 
bility. The equilibrating time is from 200000 to 400000 
MC steps/spin and the averaging time is from 500000 
to 1000000 MC steps/spin. Next, we make histograms 
at 8 different temperatures Tj(L) around the transition 
temperatures T [ f'' Tl (L) with 2 millions MC steps/spin, 
after discarding 1 millions MC steps/spin for equilibrat- 
ing. Finally, we make again histograms at 8 different 
temperatures around the new transition temperatures 
T E ' m {L) with 2 x 10 6 and 4 x 10 6 MC steps/spin for 
equilibrating and averaging time, respectively. Such an 
iteration procedure gives extremely good results for sys- 
tems studied so far. Errors shown in the following have 
been estimated using statistical errors, which are very 
small thanks to our multiple histogram procedure, and 
fitting errors given by fitting software. 

We note that only v is directly calculated from MC 
data. Exponent 7 obtained from ^ max and v suffers little 
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TABLE I: Maxima and temperatures at the maxima of V„ (n = 1, 2), C v and \ for various L with N z — 7. 



L 




max 

X. 


ymai 




Tc(C™ ax ) 


T c ( X max ) 


T c (V{ nax ) 


T c (vr ax ) 


30 


2.21115658 


25.29532589 


164.05948154 


275.71581036 


4.19027500 


4.22755000 


4.24277500 


4.24900000 


40 


2.36517434 


41.20958927 


219.30094769 


368.72462473 


4.19305000 


4.21895000 


4.23025000 


4.23500000 


50 


2.50496719 


60.82008190 


275.66203381 


463.17327477 


4.19275000 


4.21340000 


4.22210000 


4.22600000 


60 


2.59177903 


82.96529587 


329.65536262 


554.47606570 


4.19270000 


4.20940000 


4.21710000 


4.22045000 


70 


2.70129995 


109.00528127 


387.47245040 


651.24905512 


4.19250000 


4.20640000 


4.21260000 


4.21530000 


80 


2.76931676 


138.78113065 


443.00488386 


743.61068938 


4.19220000 


4.20410000 


4.20965000 


4.21205000 



errors (systematic errors and errors from v) . Other expo- 
nents are obtained by MC data and several-step fitting. 
For example, to obtain a we have to fit C™ ax of Eq. [T3] 
by choosing Co, C\ and by using the value of v. So, in 
practice, in most cases, one calculates a or f3 from MC 
data and uses the Rushbrooke scaling law to calculate 
the remaining exponent. 

Now, similar to the discussion given in subsection lIIB| 
if we work at a distance away from T C (L — oo,N z ) we 



should observe " effective critical exponents" . This is the 
case because in the finite size analysis using the multi- 
ple histogram technique, we measure the maxima of V n , 
Cy and x which occur at different temperatures for a 
finite L. These temperatures, though close to, are not 
T C (L = oo,N z ). To give a precision on this point, we 
show the values of these maxima and the corresponding 
temperatures for N z = 7 in Table |TJ For the value of 
T C (L = co, N z = 7), see Table HI 



Given this fact, we emphasize that calculations using 
Eqs. (fT2^ -(fl~7 |) will give effective critical exponents except 
of course for the case N z = 1 where the results correspond 
to real critical exponents. 

We show now the results obtained by MC simulations 
with the Hamiltonian |lj. We have tested that all expo- 
nents do not change in the finite size scaling with L > 30. 
So most of results are shown for L > 30 except for v 
where the lowest sizes L = 20, 24 can be used without 
modifying its value. 

Let us show in Fig. [T] the layer magnetizations and 
their corresponding susceptibilities of the first three lay- 
ers, in the case where J s = 1. It is interesting to note 
that the surface layer is smaller that the interior layers, 
as it has been shown theoretically by the Green's func- 
tion method a long time ago i 15 ^ 16 The surface spins have 
smaller local field due to the lack of neighbors, so thermal 
fluctuations will reduce more easily the surface magneti- 
zation with respect to the interior ones. The susceptibili- 
ties have their peaks at the same temperature, indicating 
a single transition. 

Figure [5] shows total magnetization of the film and the 
total susceptibility. This indicates clearly that there is 
only one peak as said above. 



A. Finite size scaling 

Let us show some results obtained from multiple his- 
tograms described above. Figure [3] shows the suscepti- 
bility and the first derivative Vi versus T around their 
maxima for several sizes. 

We show in Fig. 3] the maximum of the first derivative 




FIG. 1: Layer magnetizations (a) and layer susceptibilities 
(b) versus T with N z = 5 and L = 24. 

of InM with respect to (3 — (/cbT) -1 versus L in the 
In — In scale for several film thicknesses up to N z = 13. 
If we use Eq. (fl~2")) to fit these lines, i. e. without cor- 
rection to scaling, we obtain 1/v from the slopes of the 
remarkably straight lines. These values are indicated on 
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FIG. 2: Total magnetization (a) and total susceptibility (b) 
versus T with N z — 5 and L = 24. 
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FIG. 3: (a) Susceptibility and (b) Vi, as functions of T for 
several L with N z = 11, obtained by multiple histogram tech- 
nique. 



the figure. In order to see the deviation from the 2D 
exponent, we plot in Fig. [5] v as a function of thick- 
ness N z . We observe here a small but systematic de- 
viation of v from its 2D value (yw — 1) with increas- 
ing thickness. To show the precision of our method, we 
give here the results of N z = 1. For N z = 1, we have 
l/v = 1.0010 ± 0.0028 which yields v = 0.9990 ± 0.0031 
and 7/V = 1.7537 ± 0.0034 (see Figs. © and below) 
yielding 7 = 1.7520 ± 0.0062. These results are in ex- 
cellent agreement with the exact results V2D = 1 and 
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FIG. 4: Maximum of the hrst derivative of In M versus L in 
the In — In scale. The slopes are indicated on the figure. 
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FIG. 5: Effective exponent v versus N z 



72/j = 1.75. The very high precision of our method is 
thus verified in the rather modest range of the system 
sizes L = 20 — 80 used in the present work. Note that 
the result of RefUgave v = 0.96±0.05 for N z = 1 which 
is very far from the exact value. 

The deviation of v from the 2D value when N z in- 
creases is due, as discussed earlier, to the crossover to 
3D (i is not zero). Other exponents will suffer the same 
deviations as seen below. 

We show in Fig. [S] the maximum of the susceptibility 
versus L in the In — In scale for film thicknesses up to 
N z = 13. We have used only results of L > 30. Includ- 
ing L = 20 and 24 will result, unlike the case of v, in 
a decrease of 7 of about one percent for N z > 7. From 
the slopes of these straight lines, we obtain the values 
of effective 7/1/. Using the values of v obtained above, 
we deduce the values of 7 which are plotted in Fig. [7] 
as a function of thickness N z . Unlike the case of u, we 
observe here a stronger deviation of 7 from its 2D value 
(1.75) with increasing thickness. This finding is some- 
what interesting: the magnitude of the deviation from 
the 2D value may be different from one critical exponent 
to another: ~ 3% for v and ~ 8% for 7 when N z goes 
from 1 to 13. We will see below that (3 varies even more 
strongly. 

We show now in Fig. [5] the maximum of C™ ax versus 
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FIG. 6: Maximum of susceptibility versus L in the In — In 
scale. The slopes are indicated on the figure. 



FIG. 8: ln(C7 ax - C ) versus In L for N z = 1, 3, 5, 13. The 
slope gives a/v (see Eq. I13[) indicated on the figure. 
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FIG. 7: Effective exponent 7 versus N z 



Film thickness 

FIG. 9: Effective dimension of thin film obtained by using 
effective exponents, as a function of thickness. 



L for N z — 1, 3, 5, 13. Note that for each N z we had to 
look for Co, C\ and a/v which give the best fit with data 
of C™ ax . Due to the fact that there are several parame- 
ters which can induce a wrong combination of them, we 
impose that a should satisfy the condition < a < 0.11 
where the lower limit of a corresponds to the value of 2D 
case and the upper limit to the 3D case. In doing so, we 
get very good results shown in Fig. [5J From these ratios 
of a/v we deduce a for each N z . The values of a are 
shown in Table fll] for several N z . 

It is interesting to note that the effective exponents 
obtained above give rise to the effective dimension of thin 
film. This is conceptually not rigorous but this is what 
observed in experiments. Replacing the effective values 
of a obtained above in d c Q = (2 — a) /v we obtain d e g 
shown in Fig. [9] 

We note that d e Q is very close to 2. It varies from 2 
to ~ 2.061 for N z going from 1 to 13. The 2D character 
is thus dominant even with larger N z . This supports the 
idea that the finite correlation in the z direction, though 
qualitatively causing a deviation, cannot strongly alter 
the 2D critical behavior. This point is interesting be- 
cause, as said earlier, some thermodynamic properties 
may show already their 3D values at a thickness of about 
a few dozens of layers, but not the critical behavior. To 
show an example of this, let us plot in Fig. [TU] the tran- 



sition temperature at L = 00 for several N z , using Eq. 
Q~7] for each given N z . As seen, T c (oo) reaches already 
~ 4.379 at N z = 13 while its value at 3D is 4.51,2£<2S a 
rough extrapolation shows that the 3D values is attained 
for N z ~ 25 while the critical exponents at this thickness 
are far away from the 3D ones. 

Let us show the prediction of Capehart and Fisher— 
on the critical temperature as a function of N z . Defining 
the critical-point shift as 

e(N z ) = [T C (L = oo, N z ) ~ T C (3D)} /T C (3D) (22) 

they showed that 

e(N z )^-^[l + a/N z ] (23) 

where v = 0.6289 (3D value). Using T C (3L>) = 4.51, we 
fit the above formula with T C (L — 00. N z ) taken from 
Table HI we obtain a = -1.37572 and b = -1.92629. 
The MC results and the fitted curve are shown in Fig. 
[Tin Note that if we do not use the correction factor [1 + 
a/N z ], the fit is not good for small N z . The prediction 
of Capehart and Fisher is thus very well verified. 

We give here the precise values of T C (L = 00, N z ) for 
each thickness. For N z — 1, we have T C (L — 00, N z = 
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FIG. 10: Critical temperature at infinite L as a lunction of the 
film thickness. Points are MC results, continuous line is the 
prediction of Capehart and Fisher, Eq. (|23[) . The agreement 
is excellent. 
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FIG. 11: Effective exponent /3, obtained by using Eq. 1161 
versus the film thickness. 



1) = 2.2699 ±0.0005. Note that the exact value of T c (oo) 
is 2.26919 by solving the equation sinh 2 (2 J/T c ) = 1. 
Again here, the excellent agreement of our result shows 
the efficiency of the multiple histogram technique as ap- 
plied in the present paper. The values of T C (L — oo) for 
other N z are summarized in Table HU 



Calculating now M(L) at these values of T C (L = 
oo,N z ) and using Eq. [THl we obtain (3/v for each N z . 
For N z = 1, we have (3/v = 0.1268 ± 0.0022 which yields 
(3 = 0.1266 ± 0.0049 which is in excellent agreement with 
the exact result 0.125. Note that if we calculate (3 from 
a + 2/3 + 7 = 2, then (3 = (2 - 1.75198 - 0.00199)/2 = 
0.12302 ± 0.0035 which is in good agreement with the 
direct calculation within errors. We show in Fig. [IT] the 
values of (3 obtained by direct calculation using Eq. [TS] 
Note that the deviation of (3 from the 2D value when N z 
varies from 1 to 13 is due to the crossover effect discussed 
in subsection III CI It represents about 60%. Remember 
that the 3D value of (3 is 0.3258 ± 0.0044,2! 



Finally, for convenience, let us summarize our results 
in Table [EI] for N z = 1,3, ...,13. Except for N z = 1, all 
other cases are effective exponents discussed above. Due 
to the smallness of a, its value is shown with 5 decimals 
without rounding. 



B. Larger sizes and correction to scaling 

We consider here the effects of larger L and of the 
correction to scaling. For the effect of larger L, we will 
extend our size up to L — 160, for just the case N z = 3. 

The results indicate that larger L does not change the 
results shown above. Figure [T2T a) displays the maximum 
of Vi as a function of L up to 160. Using Eq. (fl"2"|) , i. e. 
without correction to scaling, we obtain \ jv = 1.009 ± 
0.001 which is to be compared to 1/v = 1.008 ± 0.002 
using L up to 80. The change is therefore insignificant 
because it is at the third decimal i. e. at the error level. 
The same is observed for j/u as shown in Fig. [T27 b): 
7/V = 1.752 ± 0.002 using L up to 160 instead of 7/1/ = 
1.751 ± 0.002 using L up to 80. 

Now, let us allow for correction to scaling, i. e. we 
use Eq.Ql]) instead of Eq. for fitting. We obtain 

the following values: 7/^ = 1.751 ± 0.002, B x = 0.05676, 
B 2 = 1.57554, w = 3.26618 if we use L = 70 to 160 (see 
Fig. [T3|) . The value of 7/1/ in the case of no scaling 
correction is 1.752 ± 0.002. Therefore, we can conclude 
that this correction is insignificant. The large value of a; 



explains the smallness of the correction. 



C. Role of boundary condition 

To close this section, let us touch upon the question: 
does the absence of PBC in the z direction cause the de- 
viation of the critical exponents? The answer is no: we 
have calculated v and 7 for N z = 5 in both cases, with 
and without PBC in the z direction. The results show 
no significant difference between the two cases as seen in 
Figs. [M] and [15] Wc have found the same thing with 
N z = 11 shown in Figs. \TE\ and [TTJ So, we conclude 
that the fixed thickness will result in the deviation of 
the critical exponents, not from the absence of the PBC. 
This is somewhat surprising since we may think, incor- 
rectly, that the PBC should mimic the infinite dimension 
so that we should obtain the 3D behavior when applying 
the PBC. As will be seen below, the 3D behavior is recov- 
ered only when the finite size scaling is applied in the z 
direction at the same time in the xy plane. To show this, 
we plot in Figs. IT51 and [T^l the results for the 3D case. 
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TABLE II: Critical exponents, effective dimension and critical temperature at infinite xy limit as obtained in this paper. 



N z 


V 


7 


Q 


P 


doff 


T C (L = oo,N z ) 


1 


0.9990 ± 0.0028 


1.7520 ± 0.0062 


0.00199 ±0.00279 


0.1266 ±0.0049 


2.0000 ± 0.0028 


2.2699 ± 0.0005 


3 


0.9922 ± 0.0019 


1.7377 ±0.0035 


0.00222 ±0.00192 


0.1452 ±0.0040 


2.0135 ±0.0019 


3.6365 ± 0.0024 


5 


0.9876 ± 0.0023 


1.7230 ± 0.0069 


0.00222 ± 0.00234 


0.1639 ±0.0051 


2.0230 ± 0.0023 


4.0234 ± 0.0028 


7 


0.9828 ± 0.0024 


1.7042 ± 0.0087 


0.00223 ± 0.00238 


0.1798 ±0.0069 


2.0328 ± 0.0024 


4.1939 ± 0.0032 


9 


0.9780 ± 0.0016 


1.6736 ± 0.0084 


0.00224 ±0.00161 


0.1904 ±0.0071 


2.0426 ± 0.0016 


4.2859 ± 0.0022 


11 


0.9733 ± 0.0025 


1.6354 ± 0.0083 


0.00224 ± 0.00256 


0.1995 ±0.0088 


2.0526 ± 0.0026 


4.3418 ± 0.0032 


13 


0.9692 ± 0.0026 


1.6122 ±0.0102 


0.00226 ± 0.00268 


0.2059 ± 0.0092 


2.0613 ±0.0027 


4.3792 ± 0.0034 
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FIG. 12: (a) V™* and (b) X max vs L up to 160 with N z = 3. 



Even with our modest sizes (up to L = N z = 21, since it 
is not our purpose to treat the 3D case here), we obtain 
v = 0.613±0.005 and 7 = 1.250±0.005 very close to their 
3D best known values v zu = 0.6302 ± 0.0001 from Ref. 
Hand v 3D = 0.6289 ±0.0008 and j 3D = 1.2390 ± 0.0025 
obtained by using 24 < L < 96 given in Ref. [53. 



IV. CONCLUDING REMARKS 

We have considered a simple system, namely the Ising 
model on a simple cubic thin film, in order to clarify the 
point whether or not there is a continuous deviation of 
the 2D exponents with varying film thickness. From re- 
sults obtained by the highly accurate multiple histogram 
technique shown above, we conclude that the critical ex- 
ponents in thin films show a continuous deviation from 
their 2D values as soon as the thickness departs from 




FIG. 13: x max vs L (a) from 20 up to 160 (b) from 70 up to 
160, for N z = 3. 
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FIG. 14: Maximum of the first derivative of InM versus L in 
the In — In scale for N z = 5 (a) without PBC in z direction 
(b) with PBC in z direction. The slopes are indicated on the 
figure. See text for comments. 
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Ln(L) 

FIG. 15: Maximum of susceptibility versus L in the In — In 
scale for N z = 5 (a) without PBC in z direction (b) with 
PBC in z direction. The points of these cases cannot be 
distinguished in the figure scale. The slopes are indicated 
on the figure. See text for comments. 




Ln(L) 



FIG. 16: Maximum of the first derivative of In M versus L in 
the In — In scale for N z = 11 (a) without PBC in z direction 
(b) with PBC in z direction. The slopes are indicated on the 
figure. See text for comments. 




Ln(L) 

FIG. 17: Maximum of susceptibility versus L in the In — hi 
scale for N z = 11 (a) without PBC in z direction (b) with 
PBC in z direction. The slopes are indicated on the figure. 
See text for comments. 
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FIG. 18: Maximum of the first derivative of InM versus L in 
the In — In scale for 3D case. 
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FIG. 19: Maximum of susceptibility versus L in the In — In 
scale for 3D case. 



1. This deviation stems from a deep physical mecha- 
nism: Capehart and Fisher™ have argued that if one 
works exactly at the critical temperature T C {L — oo, N z ) 
then the critical exponents should be those of 2D uni- 
versality class as long as the film thickness is finite. At 
T C (L = oo, N z ), the correlation in the z direction £ re- 
mains finite while those in the xy planes become infinite. 
Hence £ is irrelevant to the criticality. This yields there- 
fore the 2D behavior. However, when the system is away 
from T C (L — oo,N z ), as is the case in numerical sim- 
ulations using finite sizes, the system may have a 3D 
behavior as long as £ << N z . This should yield a devi- 
ation of 2D critical exponents. The results we obtained 
in this paper verify this picture. In addition, the predic- 
tion of Capehart and Fisher for the shift of the critical 
temperature with the film thickness is in a perfect agree- 
ment with our simulations. Note furthermore that (i) 
the deviations of the exponents from their 2D values are 
very different in magnitude: while v and a vary very lit- 
tle over the studied range of thickness, 7 and specially 
/3 suffer stronger deviations, (ii) with a fixed thickness 
N z ^ 1, the same "effective" exponents are observed, 
within errors, in simulations with and without periodic 
boundary condition in the z direction, (iii) to obtain the 
3D behavior, the finite size scaling should be applied si- 
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multaneously in the three directions, i. e. all dimensions 
should be allowed to go to infinity. If we do not apply the 
scaling in the z direction, we will not obtain 3D behav- 
ior even with a very large, but fixed, thickness and even 
with periodic boundary condition in the z direction, (iv) 
with regard to the critical behavior, thin films behave as 
systems with effective critical exponents whose values are 
those between 2D and 3D. 

To conclude, we hope that the numerical results shown 
in this paper will help experimentalists to interpret their 
data which are usually obtained at a finite distance from 



the critical point. It should be also desirable to study 
more cases to clarify the role of thickness on the behav- 
ior of very thin films, in particular the effect of the film 
thickness on the bulk first-order transition. 
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